Connectivity of Turing structures 
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It is well-known that in two dimensions Turing systems produce spots, stripes and labyrinthine patterns, and 
in three dimensions lamellar and spherical structures or their combinations are observed. We study transitions 
between these states in both two and three dimensions by first analytically deriving a control parameter and 
a scaling function for the number of clusters. Then, we apply large scale computer simulations to study the 
effect of nonlinearities on clustering, the appearance of topological defects and morphological changes in Turing 
structures. In the two-dimensional real space spotty structures we find some evidence of twin domain formation, 
of the kind seen in crystalline materials. With the help of reciprocal space analysis we find indication of other 
more general forms of order accommodation, i.e., eutactic local structures. Also a mechanism for the observed 
"connectivity transition" is proposed. 

PACS numbers: PACS numbers: 82.40.Ck, 47.54.+r, 05.45.-a 



I. INTRODUCTION 

Nature presents a fascinating diversity of patterns in plants, 
animals and other natural formations arising often from com- 
plex physico-chemical processes 1 1]. Alan Turing was the first 
to show that a simple system of coupled reaction-diffusion 
equations for two chemicals could give rise to spatial pat- 
terns due to a mechanism called diffusion-driven instabil- 
ity 101 ■ These so-called Turing patterns have since been pro- 
posed to account for pattern formation in many biological sys- 
tems, e.g. patterns on fish skin 01 Si' butterfly wings |5| and 
lady beetles to mention a few. However, the first experi- 
mental evidence of a Turing structure was reported by Castets 
et al. 17], who observed a sustained standing nonequilibrium 
chemical pattern in a single-phase open reactor with chloride- 
iodide-malonic acid (CIMA) reaction. CuiTently, there is an 
increasing interest to develop simple and plausible mathemat- 
ical models that could describe, at least qualitatively, these 
pattern formations Although there is a variety 

of models that could potentially produce similar patterns, a 
Turing system is perhaps the simplest one. 

The forms and variations of patterns generated by Turing 
systems have been studied by investigating the conditions 
for instability |12|, assuming inhomogeneous diffusion coef- 
ficients fTsll . and by introducing domain curvature fl4ll and 
growth fl5|l . In addition, symmetries in Turing systems are of 
great interest, since^ they might have biological relevance, see 
e.g. Refs. 1 13, 16, 17]. Recently, we have studied the effect of 
dimensionality by simulating three-dimensional Turing sys- 
tems, which displayed complex pattern formation 1 18]. While 
in two dimensions one obtains spots, stripes or labyrinthine 
patterns, in three dimensions complex shapes of, e.g. lamel- 
lae, spherical droplets and their combinations appear. 

Previously, studies of Turing patterns have typically con- 
centrated on reaction kinetic and stability aspects (see e.g. 119. 
I20I1 '). while the issue of pattern structure and its connectivity 
has received less attention. Here, we focus on connectivity 
of Turing patterns and its dependence on the parameters of 
the system. We present a simple way to quantitatively char- 
acterize Turing structures and their connectivity. These meth- 



ods are not only able to characterize the structures, but also 
to explain some of the effects of nonlinearities behind Turing 
patterns. We also demonstrate the existence of a "connectiv- 
ity transition" when the nonlinear interactions of the Turing 
model are varied. 

It was reported earlier that in generic 2D (3D) Turing mod- 
els nonlinear cubic interactions favor stripes (lamellae)^ with 
both morphogens being connected through the system 14,1811. 
On the other hand, quadratic interactions favor spots (spheri- 
cal structures), in which case only one of the morphogens is 
connected. Thus, a connectivity transition appears when the 
nonlinear interactions are varied from being cubic to predom- 
inantly quadratic or vice versa. In order to characterize this 
transition between separated and connected structures, we de- 
fine a dimensionless control parameter and a scaling function 
for the number of clusters. This is tested for a variety of sys- 
tem sizes and unstable modes. To characterize the resulting 
morphologies, we calculated the spatial Fourier spectrum to 
get the reciprocal space representation of the chemical con- 
centrations and followed its evolution through the connectiv- 
ity transition. 

This paper is organized as follows. Next, we briefly de- 
scribe the reaction-diffusion model of the Turing kind. Then, 
we discuss the concept of connectivity in these systems and 
the methods for characterizing the transition between differ- 
ent patterns. In Sec. lIVI we present results of comprehensive 
numerical simulations, which is followed by SecMwith a re- 
ciprocal space analysis. Then in Sec. lVII we draw conclusions. 



II. THE MODEL 

A Turing system models the evolution of the concentrations 
of two chemicals, or morphogens, and is in general repre- 
sented by the following reaction-diffusion equations 



Ut 
Vt 



DuV^U 
DvV^V 



f(U,V) 

-g{u,v), 



(1) 



where U = U{x, t) and V = V{x, t) are the morphogen con- 
centrations, and Djj and Dy the corresponding diffusion co- 
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efficients setting the time scales for diffusion. The reaction 
kinetics is described by the two nonlinear functions / and g. 

In this study we focus on the generic Turing model in- 
troduced by Barrio et al. f7], in which the reaction kinetics 
was developed by Taylor expanding the nonlinear functions 
around a stationary solution {Uc,Vc)- If terms above the third 
order are neglected, the system reads as follows 

Ut — DS\/^u + au{l — riv'^) + v{l ~ r2u) 

Vt = + v{f3 + ariuv) + u{-f + r2v), (2) 

where u ~ U~Uc and v — V—Vcare the concentration fields. 
The parameters ri and r2 set the amplitudes of the nonlinear 
cubic and quadratic terms, respectively, D is the ratio of the 
diffusion coefficients of the two chemicals, and 5 acts as a 
scaling factor fixing the size of the system. Setting D 1 
is a necessary but not a sufficient condition for the diffusion- 
driven instability to occur For details about the instability and 
the linear stability analysis of the model we refer the reader to 
Baii-io etal. HEl]. 

From the linear stability analysis PjIs*! we obtain the dis- 
persion relation and the conditions for the diffusion-driven in- 
stability as the region in fc-space with positive growth rate, 
i.e., eigenmodes u = upe^* and v — woe'*'* with eigenvalues 
A(fc) > 0. In addition, one can analytically derive the modu- 
lus of the critical wave vector 

which was here determined for the case a = —7 (set to obtain 
only one stable state at m w = in the absence of dif- 
fusion). In a discretized three-dimensional cubic system, the 
wave number is of the form 

|fc| = "^.Jiil + nl + nl, (4) 

where L is the system size and n^, Uy, are the wave number 
indices (in a two-dimensional system riz — 0). By adjusting 
the parameters and allowing only a few unstable modes, one 
can obtain several different parameter sets. As in our earlier 
work 1 18] we chose the parameters D — 0.516, a = —7 — 
0.899, [3 — —0.91 and 6^2 corresponding to a critical wave 
vector kc = 0.45, and D = 0.122, a = 0.398, (3 = -0.4 and 
6 = 2 corresponding to = 0.84. 

In this study, we vary the parameters ri and r2 in Eq. (|2|i, 
since they control the appearance of stripes or spots. By grad- 
ually changing these parameters we observe a transition from 
spotty (2D) or spherical droplet (3D) patterns to striped (2D) 
or lamellar (3D) patterns. In order to investigate this transi- 
tion, numerical simulations were carried out by discretizing 
the spatial dimensions into a square or cubic cell lattice and 
calculating the Laplacians in Eq. In all of our simulations 
we used dx = dy = dz = 1.0. The equations of motion 
were iterated in time using the Euler scheme with time step 
dt = 0.05. The boundary conditions were chosen to be peri- 
odic, and initially both chemicals were distributed randomly 
over the whole system. 



III. CONNECTIVITY 

In the numerical simulations of Eq. Q one deals with two 
concentration fields with characteristic wave lengths. In order 
to visualize this, the concentration of only one of the chem- 
icals is typically plotted with a gray scale, since in this type 
of a system the fields are in anti-phase, i.e., if there is a large 
amount of chemical U in some sub-domain, the concentra- 
tion of chemical V would be low there. These concentration 
fields vary continuously having diffuse boundaries. What do 
we mean by connectivity in patterns of chemicals? 

To answer this question one can define sub-domains domi- 
nated by either chemical U or V, provided that the amplitudes 
of the patterns is large enough. If we define the boundary 
as the interface between sub-domains dominated by different 
chemicals, we can easily locate the boundaries, since the con- 
centrations change rapidly, typically within one or two lattice 
sites. Now, if two points belong to the same domain, i.e., are 
not separated by a boundary, they are considered connected. 
The definition of the boundaries in this way is conceptual in 
the sense that in the U -dominated domains the concentration 
of V does not have to be zero, only much less than the con- 
centration of U. 

In Fig.^we show changes in the concentration fields, i.e., u 
and w, of a 2D system for different values of nonlinear parame- 
ters ri and r2 in Eq. (|2ji. When the cubic term (ri) dominates, 
the resulting stationary pattern is striped with a small num- 
ber of imperfections, see Fig. These imperfections can 
be considered as topological defects, or dislocations, which 
could serve as nucleation sites for spots. More dislocations 
appear (see Figs [l}J-C) when the strength of the quadratic 
term is made larger As the quadratic term increases, more 
spots nucleate and they arrange themselves to triangular struc- 
ture and at the same time getting rid of the remaining stripes 
(see Figs.Qf-H). Finally when the cubic term is diminished 
even further, only spots remain. In this sequence of simula- 
tions the strength of the cubic term was changed relative to the 
quadratic term by using a single control parameter P, which 
we elaborate below. Nevertheless, the transition from striped 
to spotty pattern seems to happen quite abruptly in P. In the 
last frame of this sequence (Fig.^) we see a fully stabilized 
spotty pattern with almost perfect triangular symmetry in two 
different orientations such that there is a mirror plane or twin 
boundary between them. 

Now, let us discuss patterns in Fig. ^ from the clustering 
point of view. In order to simplify this, and without loss of 
generality, we can assign zeros and ones to the whole lattice 
based on the chemical which dominates a given domain. With 
this mapping we can consider the number of clusters, which 
can be calculated using the well-known Hoshen-Kopelman al- 
gorithm i2lll as in typical percolation problems. However, 
before applying this cluster algorithm let us first visually in- 
spect possible shapes of extended clusters (Fig.Q. In Fig.^\ 
one can see that in the case of stripes the number of U- and 
^-dominated clusters is almost the same, and both types are 
extended dominantly in one of the dimensions. On the other 
hand, in the case of a spotty structure (Fig. chemical U 
appears as separate round clusters or spots, whereas chemical 
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FIG. 1 : Transition from stripes to spots. The patterns obtained after 
50000 iterations in a 100 x 100 system with kc = 0.45. Black cor- 
responds to areas dominated by chemical U (zeros) and the lighter 
color chemical V (ones). Note that the difference in parameters be- 
tween the figures is not constant: From A to I, P = 22000, 120, 75, 
65, 60, 55, 35, 15, 1. 



V forms only one connected cluster. Between these two lim- 
iting cases there is the transition region, depicted in Figs.^D- 
F, where JJ-dominated clusters appear as spots and stripes in 
the form of a "string-of-pearls". A question arises whether 
the transition from stripes to spots could be characterized in a 
reasonable way by the number of clusters? With the Hoshen- 
Kopelman algorithm we can directly determine the number of 
clusters (A^) in the system, but we also need to define a single 
control parameter that is a function of. 

One should bear in mind that many different parameter 
pairs of nonlinear terms (ri and r2) result in similar patterns. 
To obtain some general insight into the dynamics, we apply 
dimensional analysis |22] and derive a dimensionless control 
parameter for describing this behavior. The relevant dimen- 
sions are as follows 



where s denotes seconds and [c] is an arbitrary unit of concen- 
tration. Thus, the control parameter P can be written in the 
following dimensionless form 

P = n{-f. (6) 

The two limits of this parameter, i.e., P ^ oo {r2 ^ 0) and 
P ^ (ri ^ 0) yield stripes and spots (Figsn\ and ^), 
respectively, as the two extremes in 2D. The same holds for 
3D. Between these two limits a certain combination of striped 
(lamellar) and spotty (spherical droplet) patterns is expected 
to coexist (see Figs^D-Qf ). 




P 

FIG. 2: The number of U (solid line) and V (dotted line) clusters as 
a function of the control parameter P. Diamonds: kc = 0.84 in a 
100 X 100 system, squares: kc — 0.45 in a 200 x 200 system. In 
both cases the results were averaged over 20 simulations. 



IV. SIMULATION RESULTS 

In order to study the connectivity in two- and three- 
dimensional Turing patterns, we have performed extensive 
simulations using system sizes up to 5 x 10^ lattice cells and 
up to 2 X 10^ time steps to reach a stationary state. For all 
cases the results are taken as statistical averages of at least 
20 separate runs. Using P (Eq. (|6j) as a control parameter 
is plausible, since the transition from a striped (lamellar) pat- 
tern to a spotty (spherical droplet) pattern took place in a very 
naiTow region of P regardless of the individual values of the 
nonlinear parameters (ri and r2), the unstable mode (dictated 
by a and /3) or the system size. This is clearly seen in Fig.|2] 
where the number of U- and ^-dominated clusters is shown as 
a function of the control parameter P for two different cases in 
2D. In our studies the P-space was scanned by keeping either 
ri or r2 constant and by carrying out separate simulations for 
each value of P. From Fig.|2]it is clearly seen that the mode 
transition from spots to striped structure is quite sharp in P. 
The control parameter P serves the purpose of a unique tran- 
sition variable, having a well-defined value for the transition 
to occur. 

The most unstable mode for the smaller (100 x 100) system 
is with kc = 0.84, while for the larger (200 x 200) system 
it is with kc = 0.45. The larger system has smaller wave 
vector and thus the wave length, i.e., characteristic length of 
the pattern is larger. This is why the values on the vertical 
A^-axis, are of the same order of magnitude. In order to com- 
pare the numbers of clusters one can normalize it by dividing 
with A"^, where Nc — kcL/2TT, L denoting the linear sys- 
tem size (square or cube) and d the spatial dimension. A"^ 
is the maximum number of spherically symmetric clusters in 
a (i-dimensional system if the clusters were uniformly dis- 
tributed and the effect of boundaries was neglected. Due to 
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FIG. 3: The normalized number of U clusters as a function of P. 
System sizes are L — 100 (dotted), L = 175 (dashed) and L = 
500 (dash-dotted). 



FIG. 4: The normalized number of U clusters as a function of P 
in a three-dimensional system. System sizes are L = 40 (dotted), 
L = 50 (dashed), and L = 75 (dash-dot). 



the periodicity of the chemical structure, the number of clus- 
ters in the actual d-dimensional system can be estimated to 
be A^^ = (iVc + 1/2)''. However, an additional correction 
is required to take into account the effect of boundaries. One 
can estimate the number of additional partial clusters due to 
boundaries by estimating the length (area) of the boundary 
and the number of clusters within this domain (dN'^^^). As 
a result of this discussion we propose the normalization func- 
tion for the number of clusters to be 

where N{P) is the calculated number of clusters for control 
parameter P. By revising the Hoshen-Kopelman algorithm 
one could have directly calculated the number of clusters by 
taking periodicity into account, in which case Eq. reduces 
to Cd{N{P),Nc) = N{P)/N^. However, this approach was 
not implemented and the normalization was carried out by us- 
ing the Eq. 0. 

We have also studied the effect of the system size on the 
mode transition in 2D. This is shown in Fig. |3l where the 
normalized quantity C2{N{P), Nc) is plotted for J7-clusters 
against P for three different system sizes (L — 100, 175 and 
500). Neglecting the number of V^-clusters does not affect our 
conclusions, since the curves would be symmetrical as can be 
seen from Fig. |2] On the other hand, one can clearly see in 
Fig- 13 that the control parameter P succeeds in capturing the 
essential features of the transition. In addition, it can be seen 
that the normalization function of Eq. scales the number of 
clusters such that it collapses onto the same curve with only 
very small deviations outside the transition. 

If, on the other hand, one carries out the simulations for 
very small systems, finite-size effects can be observed. For 
small system sizes the C2{N{P), Nc) curve becomes very 
steep in the transition region. This would suggest that in the 



limit of small systems, the transition would become almost 
discontinuous. However, the system cannot be made infinitely 
small since the (periodic) boundary conditions start to affect 
the behavior of the system. As discussed earlier the spots tend 
to nucleate from topological defects, or dislocations, of the 
striped pattern, i.e., from the points where the stripes coincide 
(Fig.0. In the case of a small system even one dislocation can 
affect the morphology of the whole system and thus quickly 
transform stripes into a lattice of spots. In a larger system 
many dislocations have to appear at various sites to give rise 
to spots which in turn make the appearance of more spots fa- 
vorable. 

So far we have discussed our simulation results in 2D sys- 
tems. We have also studied the connectivity transition ex- 
tensively in three dimensions. In this case stripes and spots 
become lamellae and spherical droplets, respectively, and the 
structure seems more complicated especially in the transition 
region. For illustrations of three-dimensional Turing struc- 
tures we refer the reader to Ref. |18]. In 3D the transition 
does not occur at the same point with respect to P as in 2D 
since the third dimension gives to the clustering process one 
more degree of freedom, and thus it is easier for the structures 
to connect. 

This is indeed what one finds. Figure 0] depicts the nor- 
malized number of clusters for three different system sizes 
(L — 40, 50 and 75). One can see that the behavior of the 
system is different from two dimensions. Now, the transition 
occurs at a value of P which is a decade smaller than in 2D, 
since a smaller cubic nonlinear coefficient (ri) favoring lamel- 
lar structures is sufficient for increasing connectivity in three- 
dimensional space. The significance of this large drop of the 
critical P-value remains unanswered. In addition, unlike in 
2D we did not observe any finite-size effects for the smallest 
possible system sizes. 
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FIG. 5: Chemical concentration field p(r) in a system of size 200 x 
200 (top) and the corresponding diffraction patterns in fc^^fcy -space 
(bottom) for P = 22000, 150, 75, 70, 55, 12. {k^,ky) = (0, 0) is 
in the middle. 

V. RECIPROCAL SPACE ANALYSIS 

Next, we exploit another way to study the connectivity tran- 
sition, namely the discrete Fourier transform of the concentra- 
tion data in the wave vector space, k, i.e., 

p(fc)=^p(f)e*^: (8) 

This approach has been used earlier, e.g. in connection of 
reaction-diffusion systems lEsll . Turing patterns |24|, and to 
characterize the evolution of patterns |25|. The quantity p{k) 
corresponds to a diffraction pattern. Figure |5] shows a se- 
quence of the original concentration fields and their diffrac- 
tion patterns in a two-dimensional system for different control 
or transition field parameters, P. 

In Fig.|5l\, the diffraction intensity is predominantly to k^- 
direction due to more stripes in y-direction, while the smaller 
diffraction intensity around the fc^-axis is due to the stripes 
in x-direction. The distance of these diffraction peaks from 
the origin gives the length of the wave vector of the unstable 
mode, while its width in the perpendicular direction describes 
deviations of stripes from the principal directions due to stripe 
tilting and bending. In Fig. |5jJ the diffraction intensity be- 



comes more spread around the {k^, ky) diagonal, as a result 
of appearing dislocations and nucleating spots. After that in 
Figs. |5f -D the diffraction intensity starts splitting into sep- 
arate peaks due to more spots forming, then developing into 
six separate equidistant (from the origin and from each other) 
intensity peaks as evident in Figs.|5Ji-F. This hexagonal sym- 
metry in the reciprocal space is because the system evolves 
towards regular spotty pattern with predominantly triangular 
symmetry in the real space. However, these diffraction peaks 
are somewhat spread and in fact each peak split in two, which 
indicates that the triangular symmetry is not perfect over the 
whole system. 

If one examines the real space picture of spots in Fig. [Sf, 
one might think that there are fairly large disordered regions 
that deviate from the ideal triangular perfect lattice of spots. 
However, the corresponding reciprocal space picture reveals 
that this is not so, since it does not show diffuse diffraction 
pattern or a ring due to orientational disorder, but it does show 
a double peak structure as can be seen in crystalline materials 
due to twinning or low-angle grain boundaries. The existence 
of this kind of boundaries is common place in any nucleation 
mechanism and usually it is associated to the fact that twin- 
ning or low-angle grain boundaries require only small local 
displacements at low energy cost. In the present case, how- 
ever, there is no energetics involved and in Fig.ISp it is diffi- 
cult to find clear twinning, which we found in Fig.f^. There 
are other more general ways of establishing order, namely the 
appearance of eutactic local structures defined in the study 
of regular polytopes fxf\. It has been recently proposed |2^ 
that eutacticity is a very important property exhibited not only 
by crystalline and quasicrystalline lattices, but also by the ge- 
ometrical forms of some biological systems. The presence 
of order in these Turing patterns, as revealed by their clear 
diffraction patterns, suggests that the principles governing the 
preference for eutactic structures also apply to the present 
case. 

Here, we have done the reciprocal space analysis for 2D 
systems, but the same analysis can easily be extended to 3D 
systems by fixing the orientation of one of the spatial vectors, 
thus obtaining an in-plane diffraction graph. These diffraction 
patterns, however, are expected to look more or less the same 
as in 2D, which can be understood on the basis that almost 
every cross-section of a lamellae or spherical droplet pattern 
in 3D would look as striped or spotty pattern, respectively. 



VI. DISCUSSION 

In this study, we have investigated the connectivity of spa- 
tial patterns generated by the reaction-diffusion mechanism 
both in 2D and 3D. This was done by cluster analysis for the 
dominating chemical. Since several different combinations of 
parameters in a Turing system can produce similar patterns, 
we derived a dimensionless control parameter to investigate 
the effects of nonlinearities. The numerical simulations were 
consistent with the predictions drawn from the dimensional 
analysis, and the system showed a transition in the proxim- 
ity of the same P-value irrespective of the individual system 
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parameters. With help of an analytically derived normaliza- 
tion function, the number of clusters collapsed onto the same 
master curve independently of the system size or the unstable 
mode. 

From biological perspective, the characteristics of the Tur- 
ing system presented here could be important. The fact that 
the transition from striped or lamellar to spotty or spherical 
droplet structures does not depend on the size of the system 
or the specific parameters involved, but just the combination 
of the nonlinearities (P), makes the mechanism very general. 
In addition, the observation that irregular structures are not as 
stable as regular structures could be used to explain the steep- 
ness of the transition: Unfavorable structures corresponding 
to the transition domain occupy minimal volume of the phase 
space (possible structures). 

The presence of eutacticity is well known in atomic crystal 
structures. The fact that similar geometric principles seem to 
be acting also in the formation of Turing patterns could be an 
important feature in favor of their application in the study of 
the geometrical forms of some simple living organisms. This 
could also support the idea that Nature prefers geometrical 
forms generated from eutactic stars ll28il . This very simple 
analysis of Turing patterns and their morphological transitions 
due to nonlinearities could open up a new way of investigat- 
ing universal features of patterns obtained by very different 
mechanisms. 

The approach we have taken here to study pattern forma- 
tion shares common features with percolation and could be 
studied to certain extend as such. Percolation of morphogens 



to different directions is not as interesting a measure as the 
number of clusters, since it holds only a binary value and does 
not capture features of the morphological structure effectively. 
However, if one considers the field of biological applications, 
where Turing systems are often used, percolation behavior is 
interesting and has been used for studying some biological 
problems f?^. 

If one considers the patterns in Fig.[2one can observe that 
in the first frame both morphogens have percolated in the ver- 
tical direction. On the other side of the transition one of the 
morphogens has not percolated while the other has done so in 
both directions. Should one of the chemicals be favorable for 
diffusion, these two opposite stages obtained by a small vari- 
ation of nonlinear parameters could for example correspond 
to chemical concentrations in tissue as signaling via diffusion 
is either enabled or disabled. However, one should always be 
careful in making suggestions for biological applications due 
to complexities of nature not yet understood. 
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